#ifndef am_QsiAnnandale_H
#define am_QsiAnnandale_H

#include <string>
#include <math.h>
#include <tbb\concurrent_vector.h>
#include <tbb\blocked_range2d.h>
#include <boost\date_time.hpp>

#include <borealis\Raster.hpp>
#include "..\AuroraTypes.hpp"


namespace Aurora
{
	/*
	Class: am_QsiAnnandale
		Calculates the daily average incoming shortwave radiation following Annandale 2002.
		The equation uses two seperate variables depending on the location of the domain. Inland and Costal are the two choices 
		Latitude calculations are currently HARDCODED to Marmot Creek, Alberta. This WILL BE CHANGED IN THE FUTURE.

	In Variables:
		Borealis::Raster<float> Girdded DEM
		Borealis::Raster<float> Gridded Max-Min Temperature (run am_TempInterp on Dt)
		string		  Date header name
		string		  Location: "inland" or "coastal"
	Out Variables:
		Borealis::Raster<float> Qsi
	References:
		Annandale, J. G., N.Z. Jovanic, N. Benade, and R. G. Allen. 2002. Software for missing data error analysis of Penman-Monteith reference evapotranspiration. Irrigation Science 21, no. 2: 57-67.

	*/
	class am_QsiAnnandale : public AuroraTypes::ModuleBase
	{
	private:

		//in
		int DEM;
		int GRID_DT;
		int DATE;
		int LOCATION;
		int DT;

		//out
		int QSI;



	public:

		am_QsiAnnandale(std::string n)
		{
			DEM			= 0;
			GRID_DT		= 1;
			DATE		= 2;
			LOCATION	= 3;
			DT			= 4;

			QSI			= 0;

			name = n;
		}

		class CalcQsi
		{
			
			float const m_So;
			float const m_krs;
			boost::shared_ptr<Borealis::Raster<float> > const m_Dt;
			boost::shared_ptr<Borealis::Raster<float> > const m_Qsi;
			boost::shared_ptr<Borealis::Raster<float> > const m_Dem;


		public:
			CalcQsi(float So, float krs, boost::shared_ptr<Borealis::Raster<float> > Dt, boost::shared_ptr<Borealis::Raster<float> > Qsi,boost::shared_ptr<Borealis::Raster<float> > Dem)
				: m_So(So),  m_Dt(Dt), m_Qsi(Qsi), m_krs(krs), m_Dem(Dem)
			{

			}

			void operator()(const tbb::blocked_range2d<unsigned int>& r) const
			{

				
				float So =  m_So;
				float krs = m_krs;
				
				boost::shared_ptr<Borealis::Raster<float> > Dt = m_Dt;
				boost::shared_ptr<Borealis::Raster<float> > Qsi = m_Qsi;
				boost::shared_ptr<Borealis::Raster<float> > Dem = m_Dem;



				for(unsigned int x = r.rows().begin(); x!=r.rows().end();++x)
				{
					for(unsigned int y = r.cols().begin(); y!=r.cols().end(); ++y)
					{
						double dt = Dt->GetRasterAt(x,y);

						//transmittance
						double alt = Dem->GetRasterAt(x,y);
						double Tt = krs *(1+2.7*pow(10.0,-5.0)*alt) * sqrt(dt); 


						float pointQsi = So * Tt; //need to convert to MJ
						Qsi->SetRasterAt(x,y,(So * Tt)*1000000/(3600*24)); //conversion from MJ/m^2 to W/m^2
					}
				}
			}
		};

		void run(AuroraTypes::StationList* stations, AuroraTypes::VariableList* variables, AuroraTypes::InVariables& inVars, AuroraTypes::OutVariables& outVars)
		{

			//pi constant
			float PI = 3.14159F;

			//latitude in radians
			//50 N
			float phi = 0.889364355F;

			//day number
			boost::gregorian::date date = stations->at(0)->Now().GetGregorian();
			int J = date.day_of_year();

			//solar declination in radians
			//float delta = 23.5F*PI/180*sin(((2*PI)/365)*(J-80));
			float delta = 0.409*sin(2*PI/356*J-1.39);

			//118.5 MJ/(m^2 day)
			float Sprime = 118.08F;

			//max potential extraterrestrial solar radiation
			float t1 = acos(-tan(delta)*tan(phi)) * sin(phi) * sin(delta);
			float t2 = cos(phi) * cos(delta) * sin(acos(-tan(delta)*tan(phi)));
			float Drel = 1+0.033*cos(2*PI*J/365);

			float So =  (Sprime*Drel/PI)*(t1+t2 );


			boost::shared_ptr<Borealis::Raster<float> > Qsi;
			boost::shared_ptr<Borealis::Raster<float> > Dem;
			boost::shared_ptr<Borealis::Raster<float> > Dt;
			try
			{
				Dem = GetRaster(variables,inVars[DEM]);
				Dt  = GetRaster(variables,inVars[GRID_DT]);
		

				Qsi.reset(new Borealis::AsciiRaster(Dem->GetRow(),Dem->GetCol()));
				Qsi->SetXllCorner(Dem->GetXllCorner());
				Qsi->SetYllCorner(Dem->GetYllCorner());
				Qsi->SetCellSize(Dem->GetCellSize());
			}
			catch (std::runtime_error e)
			{
				throw std::runtime_error("In " + std::string(__FILE__) +" @ line " + boost::lexical_cast<std::string>(__LINE__) + e.what());
			}

			std::string location = inVars[LOCATION];
			//defaults krs to inland 
			float krs = 0.16F;
			if(location == "coastal")
				krs = 0.19F;
			

			tbb::blocked_range2d<unsigned int> range(0,Dem->GetRow(),0,Dem->GetCol());
			tbb::parallel_for(range,CalcQsi(So,krs,Dt,Qsi, Dem));

			AuroraTypes::VariableList::accessor a;
			variables->insert(a,outVars[QSI]);		
			a->second = Qsi;


		}



	};

}

#endif